Non-diffusive phase spreading of a Bose-Einstein condensate at finite temperature 
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We show that the phase of a condensate in a finite temperature gas spreads linearly in time at 
long times rather than in a diffusive way. This result is supported by classical field simulations, and 
analytical calculations which are generalized to the quantum case under the assumption of quantum 
ergodicity in the system. This super-diffusive behavior is intimately related to conservation of energy 
during the free evolution of the system and to fluctuations of energy in the prepared initial state. 
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I. INTRODUCTION 

Phase coherence is one of the fundamental properties 
of Bose-Einstein condensates. It is also a key feature in 
the present developments of the research on condensates 
which, ten years after the first experimental realization, 
go in the direction of integrating this powerful tool into 
other branches of physics, of which metrology and quan- 
tum information are two promising examples 

The problem of the condensate phase dynamics due 
to atomic interactions at zero temperature has been an- 
alyzed by different authors in theory Q and in experi- 
ment [H, 0, It is now well understood that an ini- 
tially prepared relative phase between two condensates 
will spread in time due to the corresponding uncertainty 
in the relative particle number as the relative phase and 
the relative particle number are conjugate variables. The 
phase dynamics of a two component condensate in realis- 
tic situations including harmonic traps, non stationarity 
and fluctuations in the total number of particles was an- 
alyzed in [6J, where a comparison to the experiments of 
[J] is also performed. An important conclusion was that 
the zero temperature theory could not account for the 
coherence times observed in experiment, which raises the 
question of the role of the non-condensed fraction. 

In this paper we address the fundamental problem of 
phase spreading of a Bose-Einstein condensate in a finite 
temperature atomic gas. In order to obtain simple and 
general results, we consider the ideal case of a spatially 
uniform condensate at thermodynamic equilibrium, and 
we assume that one has access to the first order tem- 
poral correlation function (aj(t)ao) of the component ao 
of the atomic field in the condensate mode. In real life, 
the situation is more complex: the atoms are trapped 
in harmonic potentials, and the measurement of phase 
coherence is a delicate procedure, usually relying on the 
interference between two condensates [4J. In the litera- 
ture two well distinct predictions exist for the long time 
spreading of the condensate phase at finite temperature, 
either a diffusive behavior (variance growing linearly in 



time) 0,1, SEl or a ballistic behavior (variance growing 
quadratically in time) [ll|. We study this problem first 
with a classical field model [H, [U, Q El Gl, E3 , where 
exact numerical simulations can be performed. We then 
explain the numerics analytically, and extend the analyt- 
ical approach to the quantum case. 

The important result that we obtain is that the vari- 
ance of the phase increases quadratically in time. This is 
at variance with the prediction of phase diffusion from the 
"quantum optics" open system approaches of @, H, H, G3 
assuming the condensate to evolve under the influence of 
Langevin short memory fluctuating forces. Our predic- 
tion results from two ingredients, (i) the system is pre- 
pared in an initial state with an energy fluctuating from 
one experimental realization to the other, here sampling 
the canonical ensemble, and (ii) the system is isolated 
in its further evolution and therefore keeps a constant 
energy. As we shall see, the combination of these two in- 
gredients prevents some temporal correlation functions to 
vanish at long times. Our prediction qualitatively agrees 
with the one of [ill ] , but not quantitatively, as we obtain 
a different expression for the long time limit of the vari- 
ance of the phase over the time squared. This difference 
is due to the fact that we take into account ergodicity in 
the system resulting from the interactions among Bogoli- 
ubov modes such as the Beliaev-Landau processes. 

In section |TT] we present the classical field model; nu- 
merical predictions for this model are presented in sec- 
tion [Till an d analytical results reproducing the numerics 
at short or long times are given in section IIV1 These an- 
alytical results are extended to the case of the quantum 
field in section [V] Wc conclude in section |VT1 



II. THE CLASSICAL FIELD MODEL 

In this section we develop a classical field model that 
has the advantage that it can be exactly simulated nu- 
merically. This will allow us to understand the physics 
governing the spreading of the condensate phase and to 
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test the validity of various approximations, paving the 
way to the quantum treatment. 

We consider a lattice model for a classical held ip{ r ) 
in three dimensions. The lattice spacings are l\, I2, I3 
along the three directions of space and dV = hhh is 
the volume of the unit cell in the lattice. We enclose 
the atomic field in a spatial box of sizes Li, L2, £3 and 
volume V = L1L2L3, with periodic boundary conditions. 
The discretized field has the following Poisson brackets 



W(ri),^*(r 2 )} 



<5ri, 



dV 



(1) 



where the Poisson brackets are such that df/dt = {/, H} 
for a time-independent functional / of the field ip. The 
field if) may be expanded over the plane waves 



(2) 



where k is restricted to the first Brillouin zone, k a £ 
[—it /la, ft/la[ where a labels the directions of space. 

We assume that, in the real physical system, the total 
number of atoms is fixed, equal to N. In the classical 
field model, this fixes the norm squared of the field: 



dVY,H(r)\ 2 =N. 



Equivalently the density of the system 

N 

P= V 



(3) 



(4) 



is fixed for each realization of the field. The evolution of 
the field is governed by the Hamiltonian 

H = J2 Eka* u a k + § £ dVf (r)f (^(^(r), (5) 



where E k is the dispersion relation of the non-interacting 
waves, and the binary interaction between particles in the 
real gas is reflected in the classical field model by a field 
self-interaction with a coupling constant g = 4irh 2 a/m, 
where a is the s-wave scattering length of two atoms. 

In general, we expect the predictions of a classical field 
model to be cut-off dependent, i.e. the predictions of our 
model may depend on the lattice spacings l a . We use 
here a refinement to the usual classical field model, which 
makes it cut-off independent for some observables like 
the condensate fraction, a quantity expected to play an 
important role here. An obvious example of a quantity 
which will remain cut-off dependent is the mean value of 
the Hamiltonian H in thermal equilibrium. 

Let us consider first the non-interacting case (<? = 0) 
in presence of a condensate. For a thermalized classical 
field the occupation numbers of the excited plane wave 
modes are given by the equipartition formula 



{a k a k ) 



E k 



(6) 



We adjust the dispersion relation Ef. in order to repro- 
duce the Bose law for the occupation numbers of the 
quantum field in the Bose-condensed regime: 



1 



2/ 37l 2 fc 2 /2m _ I 



Ek 



(7) 



For all modes with large occupation number Ek ~ 
Ti 2 k 2 /2m, while the occupation of modes with 
h 2 k 2 /2m 3> fcflT, whose quantum dynamics is not well 
approximated by the classical field model anyway, is ex- 
ponentially suppressed as in the quantum theory. 

In the interacting case, one could adapt the same trick 
of a modified dispersion relation, by including the fact 
that the relevant spectrum is not H 2 k 2 /2m but the Bo- 
goliubov spectrum (l8| . The resulting E k would now 
start growing exponentially with k when the Bogoliubov 
energy [(h 2 k 2 /2m)(2pg + h 2 k 2 /2m)} 1 / 2 reaches k B T. 

In the classical field model we restrict our analysis to 
the regime k B T ^> pg so that at energies of the order of 
k B T , the Bogoliubov energy is dominated by the kinetic 
term h 2 k 2 /2m. One can then simply use in the Hamil- 
tonian the modified dispersion relation E k as given by 
Eq.0. This is what we did in the simulations of this 
paper, so that the classical field if> evolves according to 
the non-linear equation [l!|: 



ihd t iP = {k B T exp(-/3— Aj -lj + 5 |V(r,t)| 2 | if,. 

In practice this equation is integrated numerically with 
the FFT splitting technique. 

We then introduce the density and the phase of the 
condensate mode 



«o 



(9) 



In what follows, we concentrate on three physical quan- 
tities: the condensate amplitude correlation function 

(aS(t)ao(O)), (10) 

the condensate atom number correlation function 

(SN {t) 5N (0)) where 5N = N a - (N ) , (11) 

and the variance of the condensate phase change during 
t: 

Var^(f) = (ip(t) 2 )-((p(t)) 2 where ip(t) = 0(t)-0(O). 

(12) 

The averages are taken over stochastic realizations of the 
classical field, as the initial field samples a thermal prob- 
ability distribution. 



III. CLASSICAL FIELD: NUMERICAL 
RESULTS 

We consider a gas of N = 4 x 10 5 atoms with pg = 
700 h 2 jrnV 2 ! 7, in a box of non commensurable square 
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lengths to guarantee efficient ergodicity in the system, 
in the ratio L\ ; L\ : Z| = \f2 : (1 + VE)/2 : %/3. We 
choose the number of the lattice points in a tempera- 
ture dependent way, such that the maximal Bogoliubov 
energy [(h 2 k 2 /2m)(2pg + h 2 k 2 /2m)} 1 / 2 on the lattice is 
equal to 3fcsT. 

To generate the stochastic initial values of the classi- 
cal field we proceed as follows, (i) For each realization, 
we generate a non condensed field (r) at temperature 
T in the Bogoliubov approximation as explained in [20| . 
In practice we generate complex numbers {6k} for each 
vector k on the grid according to the probability distri- 
bution 



7T Kb± 



-{\b k \ 2 e k /k B T) 



(13) 



where l k = [E k {E k + 2pg)] 1 / 2 . With a set of {6 k } for a 
given realization we build the non condensed field 



^x(r) 



/ <= ikr - p- ik r 



(14) 



where the initial value of the condensate phase 9 is ran- 
domly chosen with the uniform law in [0, 2tt[, and where 
the real amplitudes U k , V k , normalized as (7^ — V 2 — 1, 
are given by the usual Bogoliubov theory, here with the 
modified dispersion relation, so that 



U k + V k = 



E k , 



1/4 



E k + 2pg 



(15) 



(ii) We create the classical field with the constraint that 
the total number of atoms N is fixed: 



^(r) = ^L+^(r) 



(16) 



where ag = \/N — N±e ld , N± is the number of non con- 
densed atoms, 



N ± =Y,dV\iP ± {r)f 



(17) 



(iii) We let the field evolve for some time interval with 
the Eq.© to eliminate transients due to the fact that 
the Bogoliubov approximation used in the sampling does 
not produce an exactly stationary distribution. After this 
'thermalization' period we start calculating the relevant 
observables, as ip evolves with the same Eq.©. 

First we investigate the mean condensate phase change 
(ip)(t). We find a linear dependence with time, with a 
slope slightly different from the value —pg/h naively ex- 
pected, e.g. from the zero temperature Gross-Pitaevskii 
equation. The slope difference is temperature dependent 
and is expected physically to correspond to the discrep- 
ancy between the zero temperature chemical potential 
pg and the actual finite temperature one p-(T). This we 
shall confirm using Bogoliubov theory in Sec. lIVK sce also 

dH). 



In figure [TJ we show the real part of the amplitude 
correlation function of the condensate (aQ(i)ao(O)) as a 
function of time, for a temperature T = 0.17T C , where 
T c is the critical temperature of the ideal gas. The zero- 
temperature evolution e lp9t ^ h is removed so that the os- 
cillations in the figure are due to the above mentioned 
effect p(T) ^ pg. Due to the finite temperature in the 
system, the correlation function of the condensate ampli- 
tude is smeared out at long times. 
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FIG. 1: (Color online) Real part of the condensate amplitude 
correlation function (|10|) normalized to its t = value, and di- 
vided by the zero temperature evolution e !p9t//fi , as a function 
of time: in the vertical axis label, ao(i) stands for ao{t)e lp9t ^ h . 
(a) Short times behavior and (b) long times behavior. In solid 
line from an average over 500 solutions of Eq.©, in dashed 
line (red) the Bogoliubov approximation (|40p . Here the tem- 
perature is k B T = 3077.3 h 2 /mV 2/3 = 0.1711T C , where T c 
is the critical temperature k B T c = (2vr?i 2 /m)(p/C(3/2)) 2/3 of 
the ideal gas, the number of particles is N = 4 x 10 5 and the 
coupling constant is such that the zero-temperature chemical 
potential is pg = gN/V = 700ft 2 /mV 2/3 . 

Correspondingly the standard deviation of the conden- 
sate phase change increases with time, as we show in fig- 
ure for five different values of the temperature, up to 
T = 0.65T C . In all cases, at long times, we observe a 
quadratic growth of Var (p contrarily to the phase diffu- 
sion behavior cx t predicted in the literature [Tj, [8|, [9|, [10| . 

To complete the physical picture, we show in figure 
[3] the correlation function of the condensate atom num- 
ber (fTTj) . At very short times, see the beginning of the 
curves in Fig[3^,, the simulation (square symbols) con- 
firms the Bogoliubov prediction (dashed oscillating line); 
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FIG. 2: (Color online) Standard deviation of the condensate 
phase change <p(t) (|12[) as a function of time for (a) T = 
0.08245T C (lower curve) and T = 0.1711T C (upper curve), (b) 
T = 0.29467T C , (c) T = 0.453T C , (d) T = 0.6473T C . Thick 
solid line (black): numerical solution from the classical field 
model Eq.([8]) averaged over 500 realizations. Thin solid line 
(red): a linear fit. The parameters N and pg have the same 
values as in Fig[T] 



at long times, see Fig|3jD, the correlation function drops 
to a value significantly smaller than the Bogoliubov pre- 
diction (fast oscillations are not shown in the figure); a 
key point is that this long time value of the correlation 
function of the condensate atom number is not zero. 

One may fear at this stage that the classical field 
model is missing some source of damping in the dynam- 
ics of the system. However it is a well established fact 
that the classical field model is able to simulate damp- 
ing processes, including the finite temperature Beliaev- 
Landau processes [IH, [H, [H, [25|, since the interaction 
among the Bogoliubov modes is included in this model 
03 El E EE 0, El, M, M, HI HI . More quantitatively 
we now check that the damping times due to the Beliaev- 
Landau processes in the simulation are much shorter 
than the evolution times considered here. To this end, 
we extract from the simulations the temporal correlation 
functions (6*(t)6 k (0)> and (|& k | 2 (t)|6 k | 2 (0)} - (|6 k | 2 ) 2 , ob- 
tained by projecting the classical field over the corre- 
sponding Bogoliubov mode and averaging over many re- 
alizations. We show these correlation functions for the 
lowest energy Bogoliubov mode and for an excited Bo- 
goliubov mode in Fig[4j 

We come then into a paradox. On one side, the various 
Bogoliubov oscillators & k decorrelate at long times. On 
the other side, the variance of the phase change (p of 
the condensate varies quadratically at long times, which 
implies, as we shall see in Sec. IIV1 that the derivative of 
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FIG. 3: (Color online) Correlation function of the conden- 
sate atom number (JTTJ . (a) Short times, (b) long times. The 
classical field results are obtained from an average over 500 
solutions of Eq.®; they are represented by symbols in (a) 
and a solid line in (b). The dashed lines (red) are the Bo- 
goliubov approximations (|21|l (oscillating line) and ((22} (hor- 
izontal line) in (a), and only (|22p in (b). The dashed-dotted 
line (purple) is the Gaussian model. For clarity in (b) we 
washed out fast oscillations in the simulation result and in the 
Gaussian model, by averaging over consecutive points over 
a time width 0A5mV 2 ^ 3 /h. The horizontal dashed-dotted- 
dotted line (blue) in (a) and (b) is the ergodic long time limit 
prediction, described in section HVl The parameters N, T and 
pg have the same values as in Fig[T] 



the phase <fi does not decorrelate at long times, although 
it is a function of the 6 k 's; similarly, the fluctuations of 
the number of condensate atoms <5iVo, which are functions 
of the 6 k 's, do not decorrelate at long times. 



This paradox will be explained in Sec. IIVI and quanti- 
tative predictions for long times behavior of the conden- 
sate atom number correlation function and of the vari- 
ance of the condensate phase change will be derived. An- 
ticipating these analytical results, we show in Fig[5H the 
long time limit of (Var tp) 1 / 2 /t as a function of T/T c , from 
the results of the classical field simulations, but also from 
the predictions of the Bogoliubov approximation Eq. ([38|) . 
and of the ergodic theory of Sec. IIVI In figure^) we show 
the same results and predictions for the asymptotic value 
of the condensate atom number correlation function. 



5 



o 
O 

"D 
C 

CO 



o 
O 



o 
O 

"a 
c 

TO 



O 
O 




t [27tmV 2/3 /h] 



FIG. 4: (Color online) Correlation function of |fek| 2 (black 
solid line) and modulus squared of the correlation func- 
tion of 6k (green dashed line), as obtained from 2000 
realisations of the classical field simulations, (a) for the 
mode k = (0,0, 2tt/Lz) and (b) for the mode k = 
(2-k/Li, 107r/L2, 471-/L3). We define the correlation function 
of a quantity X as CorrX = (X*(t)X(0)) - \(X}\ 2 . The pur- 
ple dashed-dotted line is an exponential function of t, given 
by Eq. (124f> . It reproduces well the simulation results, and we 
have checked that the agreement is good for all vectors k ly- 
ing in an arbitrarily chosen plane in fc-space. The parameters 
N, T and pg have the same values as in Fig[T] 
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FIG. 5: (Color online) (a) Long time limit of (Var ip) 1/2 /t, 
in units of h/mV 2 ^ 3 , as a function of the temperature over 
the ideal gas critical temperature, (b) Long time limit of the 
condensate atom number correlation function, as a function of 
T/T c . In square symbols the numerical results of the classical 
field simulation (averaged over 500 realizations), in dashed 
fine (red) the Bogoliubov prediction and in dashed-dotted- 
dotted (blue) the ergodic prediction. The parameters N and 
pg have the same values as in Fig[T] 



IV. CLASSICAL FIELD: ANALYTICAL 
RESULTS 

The general procedure used here to obtain analytical 
results is the following. First one expresses the quan- 
tity of interest (the number of condensate atoms or the 
time derivative of the condensate phase) in terms on the 
amplitudes 6k of the field tp over the Bogoliubov modes, 

+V k e ^le w{ - t Hl{v,t) (18) 



where ^/>j_ is the component of tp orthogonal to the con- 
densate mode. Second one evaluates the correlation func- 
tions of products of 6k in various physical limits. 



A. Correlation function of the condensate atom 
number 

As the total number of particles is fixed, it is equivalent 
to calculate the correlation function of eSiVo in Eq. fTTj) 
and of the number of non-condensed particles N±_ . In- 
jecting the expansion Eq. (|14p for the time dependent 
non-condensed field ip± over the Bogoliubov modes into 
Eq.([T7I) we obtain 

N±(t) = lUM) + V k b*_ k (t)\ 2 . (19) 

Bogoliubov theory: In Bogoliubov theory interaction 
among the Bogoliubov modes is neglected so that at all 
times 

6k(*) = &k(0) e-*"** with uj k = £ k /fi. (20) 

As Wick's theorem applies for the initial thermal dis- 
tribution we obtain for the correlation function of the 
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condensate atom number: 
(SN (t)SN (0)) Bos = ^^[l^V^ + ^e-^ 4 ! 2 



-{U k V k f\e 



(21) 



where h k = ksT/lk is the Bogoliubov mean occupa- 
tion number of a mode for the classical field. At very 
short times, a good agreement of the Bogoliubov predic- 
tion with the simulation is observed in Fig(3^. Smearing 
out the terms oscillating rapidly at Bohr frequencies 2oj kl 
we obtain a prediction directly comparable to the coarse 
grained numerical result of Fig|3jD: 

(5N (t) 6N (0)) B o g nonosc = ^(^ 2 + V k 2 ) 2 n\ . (22) 
This amounts to considering the correlation function of 



N 



non osc / 



« = E(^ + ^ 2 ) 6 kWMi), 



(23) 



k^O 



deduced from (fT!))) by eliminating the oscillating terms 
such as &k&-k- As can be seen in Fig(3b, Bogoliubov 
theory fails at long times. Note that in the thermody- 
namic limit, where the above sum is dominated by the 
low k terms, one may approximate V k ~ —Uk, so that 
Eq.f2"2"j) is roughly half of the t = value of Eq.flHJ; in 
other words, it is approximately half of the variance of 
the condensate number. In the numerical result of Figj3l 
the correlation function drops by much more than a fac- 
tor 2. 

Gaussian theory: A possible approach to improve Bogoli- 
ubov theory consists in assuming that the 6k are Gaus- 
sian variables with a finite time correlation due to the 
Beliaev-Landau mechanism: 



\mi)h(0W=*U 



2 = f, 2 p- 2 fk\t\ 



(24) 



where 7^ is calculated with time dependent perturbation 
theory including the discrete nature of the spectrum as in 
[20| . This amounts to weighting each term of Eq. (f2"2"|) by 
exp(— 27^ \t I ) . This assumption is supported by numerical 
evidence for a single mode, see FigHJ and by an analytic 
derivation in the thermodynamic limit for one or two 
modes, see Appendix [X] Nevertheless, the resulting pre- 
diction for the correlation function of Nq, while looking 
promising at short times, see Figj3ji, is in clear disagree- 
ment with the simulation at long times, see FigJSjD. Since 
the assumption of a long time decorrelation of 6£(i) with 
&k(0) is physically reasonable, one may suspect that the 
Gaussian hypothesis is not accurate when a large number 
of modes are involved as for the correlation function of 
Nq. This is indeed the case, as we now show. 
Ergodic theory: A systematic way to calculate the long 
time limit of the correlation function is to assume that 
the non-linear dynamics generated by Eq.® is ergodic: 
at long times, the 6k(i)'s for a given realization of the 
field explore uniformly a fixed energy surface in phase 



space [30|| . In the Bogoliubov approximation for the en- 
ergy, this means that the £>k(t)'s sample the unnormalized 
probability distribution 



Poo({M) =8\ E- 



k^O 



(25) 



where the Bogoliubov energy E is fixed by the initial 
value of the field: 



E = E ^(0M0). 



(26) 



k#0 



First, for a given initial condition of the field, we calcu- 
late the expectation value of N± (t) as given by Eq. (JXSJ) 
over the ergodic distribution Eq. (|25p . which is equiva- 
lent to the temporal average of N±(t) over an infinite 
time interval. The terms of the form bb or b* b* have 
a zero mean, since the phases of the 6k 's are uniformly 
distributed over 2ir, according to Eq . (|2"5|) . To calculate 
the expectation value of the b*b terms, it is convenient to 
introduce rescaled variables 



Bk = i§ 



1/2 



6k. 



(27) 



According to Eq. (f2"5|) the real parts and the imaginary 
parts of all the are uniformly distributed over the 
unit hypersphere in a space of dimension 2A4, where 
M. = V/ dV — 1 is the number of Bogoliubov modes so 
that we obtain |i?k| 2 = 1/M where the overline stands 
for the average over the ergodic distribution (f2"5"|) . As a 
consequence the ergodic average of Nj_ is 



N 1 



E 



(28) 



k^O 



Note that this ergodic average depends on the t = value 
of the 6k 's via (f2"6f . 

Second, we average the product N± N±(0) over the 
thermal canonical distribution for the initial values 6k (0). 
This gives the long time limit of the correlation function 
of the number of condensate atoms: 



{6N (t - +oo)<Wo(0)> erg o = -Jt 

\k#0 



I + V k 2 )n k 



(29) 

This prediction is in good agreement with the simulations 
at long times, see FigJ3J) for a fixed value of the temper- 
ature, and Figj5)D as a function of temperature. Note 
that, according to Schwartz inequality, the ergodic value 
is lower than the coarse grained Bogoliubov prediction 
Eq. (f2"2"]) , as was expected physically. 

This clearly shows that the existence of infinite time 
correlations in the number of condensate atoms is a con- 
sequence of the conservation of energy during the free 
evolution of the system. 
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To understand the failure of the Gaussian model, we 
give the ergodic prediction of the long-time limit of the 
correlation function of the Bogoliubov mode occupation 
numbers n. k = |6 k | 2 , 



(n k (t -> +oo)n k '(0))er g o - (n k )(n k ') 



n k n k < 
M 



(30) 



This long-time value is non-zero, contrarily to the Gaus- 
sian model prediction. One may argue that the value 
Eq. (|30p tends to zero in the thermodynamic limit, so 
that the error in the Gaussian model looks negligible for 
a large system. However, in calculating the correlation 
function of a macroscopic quantity such as N±, a dou- 
ble sum over the Bogoliubov modes appears, so that the 
small deviations Eq. (|3"0|) from the Gaussian model pre- 
diction sum up to a macroscopic value. In other words, 
in the calculation of a given correlation function, one is 
not allowed to take the thermodynamic limit before the 
end of the calculation. 



B. Variance of the condensate phase change 

To reproduce the approach of the previous subsection 
for the phase, one should express the phase change tp(t) 
of the condensate amplitude ao as a function of the 6k's. 
It turns out that the quantity easily expressed in terms 
of the bu's is the time derivative tp. The variance of tp is 
then related to the correlation function C of tp: 



V&rtp 



ch 



dr'Ci 



(31) 



where time translational invariance in steady state im- 
poses for a classical held that C depends only on It — r'|: 



C(\t-t'\) = {<p(t)<p(t'))-{v(t)){v(t')) 



(32) 



If C(t) — > fast enough when t — > oo then ~V&rtp grows 
linearly in time. On the other hand, if C(r) has a non- 
zero limit at long times, then Var tp grows quadratically 
in time [3l| . 

To express tp in terms of the &k's, we write the equation 
of motion for ao: 



ihd = ih{a ,H} = d a *H 



(33) 



where we used d a *ip*(r) = 1/yV obtained from Eq.([2]). 
We split ip as in Eq. (fT6|) : we eliminate the condensate 
amplitude in the resulting expression for do/ao (i) by 
using |ao| 2 = N — N±, where N± is a function of the 
6 k 's, see Eq. dTOl) . and (ii) by introducing the field [32] 



A(r) 



-it) 



ip ± (r) 



(34) 



which is a function of the 6 k 's only according to Eq. (fT4f 
This leads to 



ih— = pg 
ao 



+ 



r 

A*(r)A 2 (r) 



lA(r)! 1 



YdV _ 



(35) 



The real part of the above equation gives — TiO, which is 
also —hip. 

Restricting to a weak non-condensed fraction, we drop 
the cubic terms in Eq. (|3"5|) , to obtain [33] 



hip ~ - P5 -i|^aV[A(r)+A*(r)f 



-pg- 



1L 

2V 



^(C4 + 14) 2 | 0k + 6*_: 



(36) 



It turns out that the products 6 k 6_ k generate oscillating 
terms which do not contribute to a coarse grained time 
average. It is thus useful to define 



hlpnon osc = ~ 99 - y ^ k + 



(37) 



Bogoliubov theory: By using (|20p and Wick's theorem we 
calculate the correlation function of Eq. ([36|) . By tempo- 
ral integration we obtain the variance of the condensate 
phase change 



(Var ip) Bo g 



k^O 



o sin u)ht 
t 2 + 



Qualitatively Bogoliubov theory correctly predicts a 
quadratic growth of the variance of tp at long times. As 
we show in FigJS^,, however, it is not fully quantitative: 
it does not reproduce the value of the dephasing rate ob- 
tained from the simulations. This is not surprising as in 
the full non linear theory the 6k 's interact and do not 
follow Eq.©. 

To be complete, we also give the Bogoliubov approxi- 
mation for the correlation function of the condensate am- 
plitude ao. Neglecting the fluctuations of the modulus of 
ao, one can set 



(39) 



Dropping the oscillating terms in 6k6_k and b^b* ^ in 
tp(t), which give a small contribution, we get 



(a*(t)a (0)) Bog ~ (N Q ) J] - 



k^O 



+ i£r(U k + V k ) 2 n k t 



(40) 

The resulting expression is plotted as a dashed line in 
FigfT] against the result of the simulation. 
Gaussian theory: If we add by hand a decorrelation of 
the 6k 's and assume Gaussian statistics, we get a diffu- 
sive spreading of the condensate phase change, with the 
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variance of ip 



non osc 



growing linearly at long times: 



(Var <p) 



Gauss 



(— V 



(t/ fc +V fe ) 4 n£ 



e -2 lk t _ 1 t 

\ 

Ik 



in clear contradiction with the numerical simulations. 
This prediction corresponds to a correlation function C 
vanishing at long times, whereas the correct correlation 
function has a finite limit, see Figl6l 
Ergodic theory: as in subsection IIV Al we calculate the 
long time value of the correlation function for Lp using the 
ergodic assumption. The various steps of the calculation 
are rigorously the same as in Sec. IIV Al and lead to 



C(t — > +oo) c 



\hvJ M 



(U k + V k ) 2 h k 



This prediction is in excellent agreement with the simula- 
tions: it gives the correct asymptotic value of C, see Fig® 
and from the asymptotic expression Vart^ ~ C(+oo)£ 2 
it gives the correct values of the long time limit of 
(Var ip) 1 / 2 /t, see FigEK, as a function of temperature. 



cut-off, that is on the number of Bogoliubov modes M. in 
the simulation, as is apparent on Eqs. (|29l42p . In this sec- 
tion, we therefore adapt the previous physical reasonings 
' to the quantum field case. 



A. The quantum model 

We use a straightforward generalization of the classi- 
cal field lattice model, taking here for simplicity a cubic 
lattice, as discussed in [2^, [34|, [35| . The bosonic field tp 
evolves according to the Hamiltonian 



H 



E 



2 1.2 



fifk 
~~2m 



9a 
2 



(43) 



where dk annihilates a particle of wavevector k in the first 
Brillouin zone. The dispersion relation of the wave is now 
the usual one. The total number of atoms is fixed, equal 
to N. The coupling constant go depends on the lattice 
spacing I in order to ensure a /-independent scattering 
length for the discrete delta interaction potential among 
the particles [3J, [35| : 



0.025 
0.020 
g 0.015 

D 

£ 0.010 

o 

0.005 
0.000 



10 15 20 25 30 
x [27tmV 2/3 /h] 



FIG. 6: (Color online) Correlation function C n0 n osc of the 
quantity ip non osc defined in Eq. (l37p calculated from the sim- 
ulation (square symbols; the solid line is a guide to the eye), 
or using the Gaussian theory (purple dashed-dotted line going 
to zero at long times). In dashed line (red) the Bogoliubov 
prediction. In dashed-dotted-dotted line (blue) the long time 
prediction (I42[) of the ergodic theory. The parameters T, N 
and pg have the same values as in Fig[TJ 



QUANTUM TREATMENT: ANALYTICAL 
RESULTS 



\ 9o , 

r2) = dv^> 



(44) 



Since we consider here the weakly interacting regime, we 
can restrict to a lattice spacing much larger than the 
scattering length a so that g is actually very close to 
g = Airh 2 a/m. 

To be able to use Bogoliubov theory as we did in the 
classical field reasoning, we restrict to the low temper- 
ature regime T <SC T c with a macroscopic occupation of 
the condensate mode. We thus neglect the possibility 
that the condensate is empty, which allows us to use the 
modulus-phase representation of the condensate mode: 



a 



(45) 



djdo and where 9 is a Hermitian 'phase' 



where No 

operator obeying the commutation relation 

[N ,e] = i. 



(46) 



This allows to consider the correlation of the condensate 
atom number fluctuation SNo = No — {No) but also the 
variance of the condensate phase change cp(t) = 6(t) — 
0(0), as we did for the classical field. 



So far the classical field model was very useful in re- 
vealing the physical processes governing the long time be- 
havior of the phase and atom number fluctuations in the 
condensate. However it is not a fully quantitative theory, 
as the long time limits of the correlation functions con- 
sidered here depend on the precise choice of the energy 



B. Correlation function of the condensate atom 
number 

To predict the correlation function of 5 No, we use Bo- 
goliubov theory at short times and the quantum analog 
of the ergodic theory at long times. 
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In the number conserving Bogoliubov theory [32j, |36j , 
written here for a spatially homogeneous system, one in- 
troduces the field conserving the total number of particles 



A(r) 



-w 



(47) 



where the non-condensed field ip±_ is obtained by project- 
ing out the component of the field -tip on the condensate 
mode. The field A then admits the modal expansion on 
the Bogoliubov modes 



b * u k e 



A ( r ) 2^ " l 1 

k^O 



V 



+ blV k - 



V 



(48) 



where the real amplitudes U k , V k , normalized as Ui 



1, are given by the usual Bogoliubov theory, 



u k + v k 



h z k 2 /2m 



2pg + h 2 k 2 /2m 



(49) 



Since the total number of particles is fixed to N, it is 
equivalent to consider the fluctuations of Nq or of the 
number of non-condensed atoms 



N ± = ^dFA t (r)A(r). 



(50) 



This, together with the expansion (I48p . expresses N± as 
a function of the 6 k 's. 

The equilibrium state of the system is approximated in 
the canonical ensemble by the Bogoliubov thermal den- 
sity operator 



PBo S (T) 



1 



J Bog 



*kb]b k /k B T 



(51) 



where the normalization factor Z^ 0& is the Bogoliubov 
approximation for the partition function, and where we 
have introduced the Bogoliubov spectrum 



h 2 k 2 (hfk 



2m V 2m 



1/2 



(52) 



Bogoliubov theory: In the Bogoliubov approximation for 
the time evolution, the 6 k merely accumulate a phase, at 
the frequency 0Jk = e k /h~, similarly to the classical field 
case. From Wick's theorem one then obtains 

i({JiV (t),^V (0)}) B og = Yl n k(n k + l)(U 2 k + V 2 ) 2 

+2U k 2 V k 2 cos(2w fc i) \n\ + (n k + 1) 2 [53) 

where 



n k {T) 



cxp(e k /k B T) - 1 



(54) 



is the mean occupation number of the Bogoliubov mode 
k. Note that we have considered here the so-called sym- 
metric correlation function (as {X, Y} stands for the an- 
ticommutator XY + YX of two operators) which is a real 



quantity, equal to the real part of the non-symmetrized 
correlation function. The time coarse grained version 
of the prediction ([55]) is obtained by averaging out the 
oscillating terms, which amounts to considering the cor- 
relation function of the temporally smoothed operator 
number of non-condensed particles 



k^O 



H + V 2 ) b{ 6 k + V k 2 



(55) 



Quantum Ergodic theory: Discarding from the start the 
oscillating terms in N±, as in ([53]) . we face here the 
problem of calculating the long time limit of (A(t)A(0)), 
where A is a linear function of the Bogoliubov mode oc- 
cupation numbers, 



A : 



k^o 



(56) 



As the quantum state of the system is given by the Bo- 
goliubov approximation Eq. (|51[) . we may inject a closure 
relation in the Bogoliubov Fock eigenbasis: 



1 -pY 

Z Bog ^ 



(A(t)A(0)) = — 

(»kl 

x({n k }\A(t)\{n k }), 




(57) 



where the sum is taken over all possible integer values 
of the occupation numbers, not to be confused with the 
mean occupation numbers (|54[) . 

The non-explicit piece of this expression is the matrix 
element of A(t), which may be reinterpreted as follows: 



({n k }L4(*)|{n k })=Tr [Aa(t)] 



(58) 



where the density operator er, initially a pure state in the 
Bogoliubov Fock basis, 



<K0) = |{n k })({n k }| 



(59) 



evolves during t with the full Hamiltonian H. We know 
that this evolution involves Beliaev-Landau processes 
that will spread a over the various Fock states |{n k })- 
This evolution is complex. But we need here the long 
time limit only, in which we may assume that an equilib- 
rium statistical description is possible. Since the system 
is isolated during its evolution, we take for a(t —* +oo) 
the equilibrium density operator in the microcanonical 
ensemble 371. and we calculate the expectation value of 
A with <j(t — * +oo) as we did for the classical field model. 
The calculation can be done in the thermodynamic limit. 
As shown in the Appendix [BJ one can calculate to lead- 
ing order in this limit the difference between canonical 
and microcanonical averages. 

Here the microcanonical ensemble has an energy E = 
Eq° s + Y^jto e fc n k, where E ° s is the ground state Bo- 
goliubov energy. We introduce the effective temperature 
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T e gf such that the mean energy in the canonical ensemble 
at temperature T e g is equal to E, 

= (H Bog ){T cff ) -E=J2 ek[n k (T ef! ) - n k ] (60) 

where (. . .) stands for an average in the canonical ensem- 
ble and i?Bog is the Bogoliubov Hamiltonian. Using the 
results of Appendix [Bl one gets 



**>-M><r->-T*.JS.y^ 



PWf) (61) 



where A(E) is the microcanonical average of A at energy 
E and where the apex ' stands for derivation with re- 
spect to temperature. We further use the fact that, in 
the thermodynamic limit, for typical values of the occu- 
pation numbers rik, T e e weakly deviates from the physi- 
cal temperature T. We calculate T e g by expanding ([50]) 
up to second order in T e g — T [38[ . Evaluating (j6"Tj) with 
this value of T e g, keeping terms up to the relevant order 
[HI], gives the desired result 

({n k }\A(t - +oo)|{n k }> = {A) 

{A)' , 1 / (A)' V 



^ e k (n k - n k ) 

k^O 



HBog)' 2 \(^Bog)' 
2 



(HBo S )' 



k B T 2 



(62) 



where all the canonical averages are now evaluated at the 
physical temperature T [39j. 

It remains to inject this expression into Eq. ljFTj) . The 
resulting average over n k leads to the long time value of 
the correlation function: 



<A(+oc)A(0)> - (A) 2 = 



( (A) 1 

\(HBo g ) 



Varies (63) 



where we used Wick theorem and the property dn k / dT = 
£kfik(nk + l)/kBT 2 [40|. Using Schwartz inequality, one 
can show that this long time value of the correlation func- 
tion is less than its zero time value X^7fc"-fc("fc + !)■ To 
be complete, we present an alternative derivation of our 
prediction (I64|) in the Appendix [Ul based on results ob- 
tained in [37| . We also note that the quantum ergodic 
calculation directly leads to a prediction of the long time 
limit for the correlation function of the Bogoliubov mode 
occupation numbers, see Eq. (|A17[) . 

Replacing in Eq. ([M|) the coefficients 7^ by their ex- 
pression from Eq. ([55]) , 7& = U% + V fc 2 , we obtain the long 
time value of the condensate atom number correlation 
function in the quantum ergodic theory. Note that, in 
the thermodynamic limit, this long time value scales as 
the volume V, whereas the t = value scales as V 4 / 3 . 



C. Correlation function of the time derivative of 
the condensate phase 

As in the classical field case, we first look for an expres- 
sion of the first order time derivative of the condensate 
phase operator 9 in terms of the amplitudes of the field 
A on the Bogoliubov modes. Taking as a starting point 
in Heisenberg picture 



ih—t 
dt 



[6,H], 



(65) 



we split the quantum field in a condensate part and a 
non-condensed part, 



$( r ) = -7= +$J-(r) 



(66) 



and we insert this splitting in the expression of H. Using 
the modulus-phase representation of So and the commu- 
tation relation Eq. (|4*6"|) , we obtain, using ct a = N — N±, 



dt 



9o 
V 



9o \ - 
2V ^ 



dV 



A 



.9o \ - 



2VV 



dV 



Ny + 1/2 

A>o(A> + 1) 
1 



=A + h.c. 



AU 2 + h.c. 



(67) 



The quantity (N + 1/2)/ \J N (N + 1) is actually 

1 + 0(1/Nq) so it can to a high accuracy be replaced 
by unity. Furthermore, as we did in the classical field 
model, we now keep the leading terms in A, under the 
assumption of a weak non-condensed fraction. We can 
also replace 9 by tp under the temporal derivative, since 
9(0) is time independent. We obtain [33[ 



9o 
V 



N ■ 



1 



1 



+ YdV AtA + -A 2 + -At 2 



Taking the expectation value of this expression over 
the thermal state in the Bogoliubov approximation leads 
to an expression coinciding with the value of the chemical 
potential predicted by Eq.(103) of [35j, which includes in 
a systematic way the first correction to the pure conden- 
sate prediction pgo [41] : 



k^O 



(69) 

At this order of the expansion, this analytically shows 
that —h(d9/dt) is the chemical potential of the system. 
We now turn to various predictions for the symmetrized 
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correlation function of d(p/dt, 

*w-5<{(i#)w- (!*)»»} >-<!*>'■ <™» 

Bogoliubov theory: At a time short enough for the inter- 
actions between the Bogoliubov modes to remain negli- 
gible, one can apply Bogoliubov theory to get 



C s B ° S (r) 



(U k + V k ) 4 {n k (n k + 1) 



- (#) E 



- cos(2w fc t) [nj + (n* + l) 2 



(71) 



The temporal coarse grained version of this correlation 
function is obtained by averaging out the cosine terms, 
which amounts to considering a temporal derivative of tp 
freed from the oscillating terms bb and b'b': 



dt 



9o 



N-U 



v k ) 



k^O 



(72) 



Quantum ergodic theory: We directly apply to the 
smoothed temporal derivative (|72[) the reasoning per- 
formed in the previous subsection. Up to an addi- 
tive constant, Eq. ([72|) is indeed of the form (|56|) . with 
Ik = —{go/KV){U k -\-Vk) 2 - From [[S3]) we therefore obtain 
the long time behavior of the phase derivative correlation 
function 



J2k^o( U k + V k ) 2 e k n k (n k + 1) 



Emo 4 n k(nk + 1) 



The long time limit of the variance of the phase difference 
is then 1421 



Var ip • 



C° rgo (+oc)t 2 



(74) 



Although our conclusion of a ballistic behavior for the 
phase agrees qualitatively with [ll| , the explicit expres- 
sion of the coefficient of t 2 differs from the one of [llj due 
the fact that we account for interactions among Bogoli- 
ubov modes such as the Beliaev-Landau processes leading 
to ergodicity in the system, while in jll] the many-body 
Hamiltonian is replaced by the Bogoliubov Hamiltonian 
in the last stage of the calculation. As can be seen from 
(|73|) using Schwartz inequality, ergodicity results in a re- 
duction of phase fluctuations with respect to the Bogoli- 
ubov prediction. 

In the thermodynamic limit, analytical expressions can 
be obtained for this ergodic prediction. In the low tem- 
perature limit kgT -C pg, 



^ergo/ . \ 8tt 4 a 2 £ / k B T 



pg 



(75) 




k R T/pg 



FIG. 7: (Color online) In the regime T <§C T c for a weakly 
interacting Bose gas, quantum ergodic prediction (|73p for the 
long time limit of (Var (p) 1 ^ 2 /t, in the thermodynamic limit. 
When expressed in units of (a 2 £/V) 1,/2 pg/h, (Var ip) 1 ^ 2 /t is 
a function of ksT j pg only, that is readily calculated numeri- 
cally (solid line) or that may be approximated by asymptotic 
equivalents (|75I76|I in the low temperature or high tempera- 
ture limit (green dashed line). Note that the dimensionless 
quantity a 2 £/V may also be written as */ pa? / (N \J~^n) . 



where £ is the healing length such that h 2 /rn£ 2 = pg. 
This tends to zero at zero temperature [H]]. In the high 
temperature limit k B T 3> pg, 



pcrgo 



(+oo) 



12C(3/2) 2 a 2 A (k B T 
5C(5/2) V 



(76) 



where the thermal de Broglie wavelength obeys A 2 = 
27r7i 2 /rnksT and where £ is the Riemann Zeta function. 
Here we have identified go to g [43J. In Fig [7] we give 
the quantum ergodic prediction for limt_ (00 (Var (p) 1 ^ 2 jt 
calculated numerically, which is a universal function of 
kgT/ pg when expressed in the right units. 



VI. CONCLUSION 

We have investigated theoretically the phase spreading 
of a finite temperature weakly interacting condensate. 
The gas is assumed to be prepared at thermal equilib- 
rium in the canonical ensemble, and then to freely evolve 
as an isolated system. After average over many realiza- 
tions of the system, we find in classical field simulations 
that the variance Var ip of the condensate phase change 
grows quadratically in time. This non-diffusive behav- 
ior is quantitatively explained by an ergodic theory for 
the Bogoliubov modes, the key point being that conser- 
vation of energy during the free evolution prevents some 
correlation functions of the field from vanishing at long 
times. We have extended the analytical treatment to the 
quantum field case and we have determined the coeffi- 
cient of the t 2 term in the long time behavior of Var<^, 
see Eg. (1731 . This analytical result holds at low tem- 
perature T <C T c and in the weakly interacting regime 
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pa 3 <C 1, for a large number of thermally populated 
Bogoliubov modes, and relies on the assumption that 
the (although weak) interaction among the Bogoliubov 
modes efficiently mixes them (quantum ergodic regime). 

A physical insight in our result is obtained from the 
following rewriting 



Vax<p(£) ~ — 



t 2 ( d T p 



h 2 \dr(H) 



VariJ 



(77) 



where Var H is the variance of the energy of the gas, here 
in the Bogoliubov approximation and in the canonical 
ensemble, n(T) is the chemical potential of the system as 
given by Eq. (p5)l and (H)(T) is its mean energy in the 
Bogoliubov approximation. 

This formula may also be obtained from the following 
reasoning. For a given realization of the system, of energy 
E, the long time limit of the condensate phase change 
<p(t) can be shown to behave as 



<p.(t) ~ -p micro (E)t/h, 



(78) 



where /imicm is the chemical potential calculated in the 
microcanonical ensemble [45| . For a large system, canon- 
ical energy fluctuations around the mean energy (H ) (T) 
are weak in relative value so that one may expand 
Mmicro(-E') to first order in E — (H). Taking the vari- 
ance of <p(t) over the canonical fluctuations of E then 
leads to ([771) . since dTn/dT{H) ~ dEfimicmdH)) for a 
large system. 

This reasoning shows that a necessary condition for 
the observation of an intrinsic diffusive spreading of the 
condensate phase change is a strong suppression of the 
energy fluctuations of the gas. To this end one may try 
to prepare the system in a clever way, starting with a 
pure condensate and giving to the system a well defined 
amount of energy, e.g . by a reproducible change of the 
trapping potential [461 ] . Alternatively one may try to fol- 
low a given experimental realization of the system, mea- 
suring the phase of the condensate in a non-destructive 
way and replacing ensemble average by time average. 
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APPENDIX A: TEMPORAL CORRELATION 
FUNCTION OF THE BOGOLIUBOV MODE 
OCCUPATION NUMBERS 

Using the master equation approach developed in 
quantum optics [13, EH , we calculate the temporal cor- 
relation function of the operator 6 q 6 q giving the number 



of Bogoliubov excitations in the mode of wavevector q, 
in the thermodynamic limit and including the Bcliacv- 
Landau coupling among the Bogoliubov modes. This is 
useful to motivate the Gaussian model introduced in sec- 
tion IIV1 and to estimate the time required for the cor- 
relation function (A(t)A(0)) — (A) 2 , where A is of the 
form Eq. (|56p . to depart from its value predicted by the 
Bogoliubov theory. 

The idea of the master equation approach is to split 
the whole system in a small system S and a large reser- 
voir R with a continuous energy spectrum. Treating the 
coupling W between S and R in the Born-Markov ap- 
proximation one obtains a master equation for the den- 
sity operator a$ of the small system. Here the small 
system is the considered Bogoliubov mode, with unper- 
turbed Hamiltonian H$ = eq6^6 q , and the reservoir is 
the set of all other Bogoliubov modes, with unperturbed 
Hamiltonian Hr — H^og — H$. In the thermodynamic 
limit, the reservoir indeed has a continuous spectrum, 
whereas the small system has a discrete spectrum. The 
coupling W between S and R is obtained from the next 
order Bogoliubov expansion of the Hamiltonian, that is 
from the part of the Hamiltonian cubic in the field A, 



9oP 



1/2 



J2dVA\A + At)A. 



(Al) 



Inserting the modal decomposition Eq. (l48[) in i? C ub, we 
isolate the terms that are linear in 6 q [49(: 



W = 99 



1/2 [S q flt 



Vlr 



(A2) 



k,kVO,q 

+ 5 q ,k+k'Sfe,fc'£k&k' + <5 q) k/-k2Cfc ) k/& k &k' (A3) 

where the operator R acts on the reservoir only and the 
coefficients have the explicit expressions: 

A k , k , = U q V k V k , + (U q + V q ){U k V k , + U k ,V k ) + V q U k U k , 
B k , k . = U q U k U k , + (U q + V q ){V k U k , + U k V k .) + V q V k V k , 
c k , k > = u q v k u k , + (U q + V q ){U k U k , + V k V k .) + v q u k v k ,. 

As a consequence of momentum conservation for the 
whole system, the action of R (respectively R') changes 
the reservoir momentum by —fi(\ (respectively ?iq). 

Let us denote with a tilde the operators in the interac- 
tion picture with respect to the Bogoliubov Hamiltonian 
Hs + Hr. In the Born-Markov approximation [13, HH the 
master equation, for the density operator of the small sys- 
tem in contact with the reservoir in an equilibrium state, 
reads 1501 



^Tr R {[W(t),[W(t- 



where Tr^ denotes the trace over the modes of the reser- 
voir and the equilibrium density operator of the reser- 
voir <r^ q is supposed here to be the Bogoliubov thermal 
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equilibrium at temperature T . We expand the double 
commutator; because of momentum conservation, the re- 
sulting terms that contain two factors 6 q or two factors 
6 q exactly vanish when one performs the corresponding 
traces over the reservoir. Coming back to Schrodinger's 
picture we finally obtain 



= ■4[e g &q6q,0- S f]+r 9 SqC7sS^H-r+S^(7sf6, 



^{r-S q Sq + r+SqS q ,^} 



q 

(A5) 



where { , } is the anticommutator and the new mode fre- 
quency is l q — e q + hA q . The effect of the reservoir on the 
small system is then characterized by a frequency shift 
A q of the mode, whose explicit expression we shall not 
need here (5lj . and by two transition rates and T~ 
given by the Fourier transform of reservoir correlation 
functions at the mode frequency: 

2 r+oo 

K = dre-^^ n TT R [R\T)Ra^} (A6) 

il J — oo 
2 r+oo 

rY = -TT / dre^^ n Tr R [R(T)R^a^}. (A7) 



Since the reservoir is here at thermal equilibrium, the two 
rates are not independent but T~ = 
from the Bose law property 1 + 
are then conveniently characterized by their difference 
r, = r- - T+. One finds 



This results 
Ufe. The rates 



r _ 9%P 
1 q — 



d 3 k [4Cl k ,(n k - n k ,)S(e q + e k - e k ,) 



+2B 2 Kk ,{l + n k + n k ,)6(e k + e k , - e g )] (A8) 

where k' stands for |k — q| in the integrand. From [24| 
one checks that T q is simply the standard Beliaev-Landau 
damping rate for the Bogoliubov mode q, the contribu- 
tion in C 2 corresponding to the Landau mechanism and 
the one in B 2 to the Beliaev mechanism. 

We now proceed with the calculation of the temporal 
correlation function of two operators As, Bs of the small 
system, the whole system being at thermal equilibrium. 
The quantum regression theorem [53[ states that 

(A s {t)B s ) = ((A s ))(t) ^ Tr s [A s af(t)} (A9) 

for t > 0, where the effective density operator a s s is in 
general not hermitian nor of unit trace but evolves with 
the same master equation as as with the initial condition 



of (0) = Bsvl 



(A10) 



where a c ^ = e~@ Hs /Zs is the unit trace equilibrium solu- 
tion of Eq. (|A5[) . Using the invariance of the trace under 
a cyclic permutation we obtain 



dt 



((As)) 



-^(([blAs^ + bliAsM)) 
S (( [S q ,A s ]S q + 6 q [ J 4 5 ,6 q])> .(All) 



Specializing to As = B' s = fe q or S q and As = Bs — 
riq = 6 q 6 q leads to linear first order differential equations 



for ({As))(t) that are readily solved: 



(& q (^ q > 



,(a q -r q /2)t 



{ftq + i) e (-^-iV2)t ( A13 ) 



(nq(t)n q ) - n 2 = n q (n q + l)e 



-r a t 



(A12) 
(A13) 
(A14) 



In the classical field limit, where n k + 1 is assimilated 
to n k , this justifies the Gaussian theory of section [TV] 
In both the classical and quantum cases, this shows that 
the occupation numbers decorrelate with the rates T q 
corresponding to the Beliaev-Landau processes. These 
rates have a non-zero value in the thermodynamic limit. 

The present calculation is readily extended to the in- 
clusion of two Bogoliubov modes in the small system, of 
wavevectors q and q' ^ q. The coupling of the small 
system to the reservoir now takes the form 



W 2 = gp 



1/2 



6 q i? q 



h.c. 



(A15) 



where the operators R have the same structure as in the 
single mode case, except that the double sum over k, k' is 
restricted to values different from 0, q, q'. In the result- 
ing master equation for the density operator of the two 
modes, the only issue is to see if there will be crossed 
terms between the two modes, involving e.g. the product 
of & q with bqi. By calculating the trace over the reser- 
voir of the corresponding product of operators R, e.g. 
Trij[i? q (r)7? q ,f7^ q ], we find in general that all crossed 
terms vanish, because of momentum conservation [52j . 
The master equation therefore does not couple the two 
modes, and one obtains 



,(*) 



n„u n„>) - n, 



,7V = 0, for q' ^ q, (A16) 



as is assumed in the Gaussian model for the classical field 
of section ITVl 

It is instructive to compare the long time limit of the 
predictions Eqs. (|A14| IA16j) to the quantum ergodic pre- 
diction. Adapting the reasoning leading to Eq. fMl) . we 
obtain the quantum ergodic result 



(n q (+oo)n q /) - n q n q 



e q e q >n q {n q + l)n q ,{n q > + 1) 



(A17) 

In the thermodynamic limit this tends to zero, as in the 
master equation approach. 
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APPENDIX B: DEVIATION OF 
MICROCANONICAL AND CANONICAL 
AVERAGES 

We wish to calculate the thermal expectation value of 
an observable A in the microcanonical ensemble rather 
than in the canonical one. For convenience, we shall 
parametrize the problem by the temperature T of the 
canonical ensemble. Restricting to the thermodynamic 
limit, where k B T is much larger than the typical level 
spacing of the system, we calculate the first order devia- 
tion of the two ensembles. 

We start with the usual integral representation of the 
canonical ensemble in terms of the microcanonical one: 



(A)(T) = 



JdEA{E) 



,S(E)/k Be -0E 



-f)E 



(Bl) 



where the density of states is written in terms of the 
exponential of the microcanonical entropy S(E), A(E) 
and (A) (T) stand for the expectation value of A in the 
microcanonical ensemble of energy E and in the canonical 
ensemble of temperature T respectively, and (3 — l/k B T. 

In the thermodynamic limit we expect the integrand 
to be strongly peaked around the value Eq(T) such that 



d_ 

dE 



S(E) 

k B 



f3E 



S'[E (T)]/k B - 13 = 0, 



E=E (T) 



(B2) 

where f'(x) stands for the derivative of a function / with 
respect to its argument x. We then expand u(E) = 
S(E)/(k B ) — f3E up to third order in E — E a and we 
approximate the integrand as 



ME) 



e u(E Q ) e S" (E a )(E~E ) 2 /2k E 



1 



(E-EofS^iEo)/^ 



.(B3) 



We also expand A(E) up to second order in E — Eq. 
Performing the resulting Gaussian integrals leads to 



(A)(T)-A[E (T)] = 

s (3 He ) , 



kt 



A'(E ) 



\S"(Eo)\ 



2\S»(E )\ 

A" {Eq) 



(B4) 



This relation can be inverted to first order, to give the mi- 
crocanonical average as a function of the canonical one; 
to this order, we can assume that A[Eq(T)] = (A)(T) in 
the right hand side of (|B4j) . Furthermore, using the im- 
plicit equation (|B2|) one is able to express the derivatives 
with respect to E in terms of derivatives with respect to 
T, e.g. S"[E {T)} = -1/[T 2 £^(T)]. This leads to 



A[E (T)]-(A)(T) = -k B T 



(A)'(T) T{A)"(T) 



E' {T) 2E' Q (T) 



(B5) 

It is actually more convenient to parametrize the result 
in terms of the mean canonical energy (H ) (T) rather than 



in terms of E (T). Applying (|B5[) to A = H allows to 
calculate Eq(T) — (H)(T) to first order. One then uses 
the first order expansion 



A[(H)(T)] 



A[E (T)] - 
1 d 

Em dr 



^[(H)(T)-E (T)]x 
{A[E (T)]}+... (B6) 



In the first order term of this expression, we can replace 
A{Eq(T)] by the canonical average (A)(T), and we can 
identify Eq(T) with (H)(T); we can do the same identi- 
fication in the right hand side of (|B5|) . We obtain [54[ 



A[(H)(T)]-(A)(T) 



1 

—k B T' 
2 



d fd{A)/dT\ 
~dT \d(H)/dTj 



(B7) 



APPENDIX C: ALTERNATIVE DERIVATION OF 
THE LONG TIME LIMIT OF CORRELATION 
FUNCTIONS 

We present in this section an alternative derivation of 
the ergodic result for the correlation function of an 
hermitian operator A, here introduced in (|56p . The long 
time limit of the correlation function is rigorously defined 
in terms of the temporal average 

C A (+^) = lim - f dr [(A(r)A(O)) - (A) 2 } . (CI) 



We then insert in (|C1[) a closure relation using the exact 
iV-body eigenstates |m) of the interacting system with 
eigenenergies E m . In the absence of degeneracies we ob- 
tain a single sum over m, 



C A {- 



•00 



E 



p m {m\A\mf 



^ Pm (m\A\m) 



(C2) 

Here the p m = Z 1 exp(— f3E m ) are the statistical weights 
defining the average in the canonical ensemble. Equation 
(|C2j) . specialized for jk = {go/V){Uk + Vk) 2 , is equivalent 
to Eq.(22) in [llj for the dephasing time, provided one 
replaces there H' by A. This makes the link between our 
approach and the one of [ll| . 

The delicate point is now to relate the formal ex- 
pression (|C2j) (involving the unknown exact eigenstates 
|m)) to an explicit expression treatable in the Bogoli- 
ubov approximation. If one directly approximates the 
exact eigenstates by eigenstates of the Bogoliubov Hamil- 
tonian, \m) ~ |{nk}), as done in [ll[ (see Eq.(61) there), 
one obtains the Bogoliubov result 



-oo) 



]ltn k {n k + 1), 



(C3) 



which is a good approximation for the t — value of the 
correlation function, but not for its long time limit. We 
argue that the exact eigenstates are in fact coherently 
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spread over a large number of Bogoliubov eigenstates of 
very close energies, because of the Beliaev-Landau cou- 
plings among them. Following [33], we thus assume that 

(m\A\m) ~ A(E m ) (C4) 

where A(E m ) is the microcanonical ensemble average at 
the energy E m , a thermodynamic quantity that is now 
treatable in the Bogoliubov approximation as we have al- 
ready done in Eq. (frj2|) [56| . After average over the canon- 



ical distribution for the energy E m = E, we then obtain 
for the correlation function, 



C A (+oo)-([A(E)-(A)] 2 )c,(^^y S j VartfBog, 

(C5) 

where (...) stands for the canonical average at tempera- 
ture T. We recover Eg. ([55)) . 
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T ■ 



£ 



n k (T)} 



(^Bog)'(T) 
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